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Abstract 

We quantify the correlation between earthquakes and use the same to distinguish between rele- 
vant causally connected earthquakes. Our correlation metric is a variation on the one introduced 
by Baiesi and Paczuski (2004). A network of earthquakes is constructed, which is time ordered and 
with links between the more correlated ones. Data pertaining to the California region has been 
used in the study. Recurrences to earthquakes are identified employing correlation thresholds to 
demarcate the most meaningful ones in each cluster. The distribution of recurrence lengths and 
recurrence times are analyzed subsequently to extract information about the complex dynamics. 
We find that the unimodal feature of recurrence lengths helps to associate typical rupture lengths 
with different magnitude earthquakes. The out-degree of the network shows a hub structure rooted 
on the large magnitude earthquakes. In-degree distribution is seen to be dependent on the density 
of events in the neighborhood. Power laws are also obtained with recurrence time distribution 
agreeing with the Omori law. 

Key Words: Network Analysis, Earthquake Networks, Earthquake Correlations, Earthquake 
Recurrences, Nonlinear Dynamics. 
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I. INTRODUCTION 



New ways to analyze earthquake catalogs have emerged in recent times in the wake 
of an objective approach where the subjective/arbitrary element in the identification of 
aftershocks, foreshocks etc. have been discarded g], [lO, H, Q. For example, attempts 
to identify aftershocks, foreshocks etc. using correlations between events based entirely on 

and aftershock characterizations, on the same footing is another approach [2|, |12|, |15| which 
reveals some basic characteristcs that affect them all, and which depends only on magnitude 
and, perhaps, on finite size effects. That the aftershocks and main shocks are not caused by 
different relaxation mechanisms and no clear distinction can be made between them is getting 
accepted by a significant number of researchers jlO . 
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2l[ |. All this leads to identification 



of interesting new facets of earthquake dynamics and patterns arising therefrom. 

The clustering of earthquakes in space and time suggests that events that follow in time 
are, to a certain extent, causally related to the earlier ones. However, restricting the causal- 
ity connection to the somewhat arbitrary mainshock-aftershock scenario may not be enough. 
Rather, the causality connection can be extended to a cluster of events that are strongly cor- 
related based on data analysis. We adopt this approach herein and study spatial recurrences 
by constructing a network of events that are correlated above a certain threshold. While a 
unique correlation relation may not be possible, we follow Baiesi and Paczuski [lO| and base 
our correlation relation on the two main statistical laws that characterize earthquake data — 
the Gutenberg- Richter law lisf and the fractal distribution of earthquake epicenters (see, 
for example, 
> m vary as 



31|). The former law states that the number, A^, of earthquakes of magnitude 

N{m) ~ 10-"™ (1) 



where 6 is a constant ~ 1, but does vary a little with the region and the catalog (see, for 
example, ISj). The latter is based on fractal analysis and is of recent origin, but has been 
shown to be quite robust through the analysis of various datasets from various regions of 



the globe (see, for example, |l9l. l3l|): 



A^'(/) oc r 



(2) 
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where A^' is the number of pairs of points separated by a distance / and df is the fractal 
dimension. 

If we combine the above two laws, we may state that the average number of events that 
can occur in the region within a distance lij separating two events i and j is 

= Kl -1 10'^""^ Arm (3) 

where K is a. constant of proportionality that may be related to the seismic activity of the 



region {q], and Arrii is the accuracy in measurement of the magnitude mj. Note that the 
events of the catalog are assumed to be time ordered with j > i and tj > ti. We may then 
define a correlation relation between earthquakes as 

= (4) 

and note that the correlation between two events is a maximum when riij is at a minimum. 

We construct a growing network by connecting the correlated events together. For this 
purpose, correlation thresholds, Cth, are defined and only those event pairs with Cjj > Cth 
are linked together in the network. The basic premise is that all subsequent events in the 
catalog are recurrences to the selected event, but a pruning is being done based on Cth to 
ensure that only the most correlated are selected. 

Our metric for measuring correlation is similar to the one used by Baiesi and Paczuski 



10| (hereafter referred to as BAPA), but, since we are focussing on recurrences rather than 
aftershocks, we have chosen to omit the time factor tij (= tj — ti) from the expression 
occurring in BAPA. The metric used here ensures that spatially closeby events are more 
correlated and, also, a higher magnitude event contributes a larger value to correlation than 
a lower magnitude one at that same location. In particular, our metric does not, unlike in 
BAPA, leave out events that are separated by large time intervals, as long as they satisfy 
the correlation relation, Eq. |H 

It is to be also noted that, in BAPA, the question that was asked was as to which previous 
event the present event is an aftershock of. To this end, Eq. [3] (with an added factor of tij 
in the equation) was evaluated in BAPA by holding j fixed and varying i from the set of 
all previous events. In such a procedure, rrii as well as I'^j keep changing in each evaluation. 
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In our case, we go forwards in time and ask which event is returning to the same location. 
Therefore, i and rrii are fixed, j varies between evaluations and, only lij keeps varying. 
Thus, in our case, each event selects a spatial window proportional to its magnitude and the 
subsequent events closer to it get more weightage in the computed correlation value. In the 
case of BAPA, the size of the spatial window depends on, not the magnitude of the event 
whose associations are being looked for but, the earlier events which are being sampled to 
identify associations. Our procedure identifies recurrences whereas aftershocks are identified 
in BAPA. 



A. The Catalog 

Our results have been based on a study of a seismic event data catalog of California region 
(for the rectangular region covering (120.5°W-115.0°W) longitude and (32.5°N-36.0°N) lat- 



itude), for the period January, 
to this region |2|-y, U, U lul, \l2 
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984-Deceniber, 2002 [l| . There are many studies pertaining 
, Q and this makes it possible to critically 



compare results. The catalog is assumed to be homogeneous for this period and is complete 



for events with m > 2.5 



33| . The relative location errors are of the order of 100 m |27| . 



The number of events in the catalog is 19569. The b value for the catalog obtained from a 
linear fit to the plot of logA^ vs. m of Eq. [T] is 0.97; the linear fit was good in the range 
2.5 < m < 6. 

Sub-catalogs with higher magnitude thresholds such as mth = 3, 3.5 and 4.0 and subsets 
for shorter periods were also employed to check the dependence on m and on the number of 
earthquakes. 

II. METHODOLOGY 



The absolute value of K is not critical for this study since only the relative correlation 
values matter in separating the most correlated from the less correlated. We have chosen 
to set K = 10~^ as was done in BAPA. Again, we have taken df = 1.6 as in BAPA which 
is a good representative value given the many estimates available 0, Am, = 0.1 as 
available from the data site [l| . 

The distribution of correlation values between all i-j pairs in the data set shows a scale 
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FIG. 1: The log-log plot of the probability distribution, P{c), of the correlation values (c) between 
seismic events, versus c. The latter has been binned logarithmically and are between pairs of 
events in the network satisfying m > =2.5. The linear fit in this log-log plot shows that the 
distribution is a power law with an exponent of ~ —1.7. A power law distribution is observed for 
the other rrifh values too. The range in c values decreases with increase in rrith- 



u 
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free distribution. This is evident in the log-log plot of -P(c) vs c (Fig. [T]), where c is the set 
of binned correlation values and -P(c) the probability values of those bins. The linear regime 
in this plot ranges over 8 orders of magnitude. The slope of the linear region is ~ —1.7 
for the catalog subset with rrith = 2.5, and (not shown here) —1.6 for rrith = 4.0. The 
broad distribution implies that, typically, an event i has a few events j which are strongly 
correlated to it and a larger number of j's which are weakly correlated to it. By specifying 
correlation thresholds Cth, we may select only those i-j pairs from the data set that have Cij 
> Cth- Outgoing links are drawn from each event to another if > Cth- Clusters result, 
where each cluster consists of an event and its highly correlated recurrences. A cluster can 
have as few as two connected events. There may also be isolated events with no clusters 
around them. By varying Cth, networks of different levels of density are generated. 

Since many of the commonly observed patterns in seismic events data have magnitude 
dependence, networks were also generated by selecting only those events with magnitudes 
greater than a threshold value, rrith- In this case too, correlations are calculated between 
the events of the subset of the data points with m > rrith and, subsequently, an appropriate 
Cth imposed for generating networks with links between the qualifying pairs of events. 

We have operated in the range of 10'' < Cth < 10^°. Higher thresholds result in too 
sparse a network. While Cj,- < 10^ was also obtained with many i-j pairs in the data 
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FIG. 2: Network constructed by three different procedures are presented for comparison sake. 
Fig. [2]Ja) has been constructed using the procedure advised in DAGPA, Fig. Wii>) has been con- 
structed using the procedure of BAPA and, finally, Fig. Wic) employs our own procedure. Same 
data set has been employed for all the three panels and the epicenter of Hector Mine earthquake has 
been selected as the originating node of aftershocks/recurrences. Procedure of BAPA is restrictive 
in time while the procedure of DAGPA is minimalistic and selcts an event only if it is closer than 
the previously selected event. Our procedure selects the full set of recurrences with aj > cth- 
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set, we have verified that lower Cth values do not yield any new information while at once 
making the network too dense and unwieldy. In Fig. [21 we present a part of the network 
for visualization purposes. Three panels have been provided for ready comparison of the 
differences in the networks generated by three different procedures, viz. by a procedure 
developed by Davidsen, Grassberger and Paczuski [l^i] (hereafter referred to as DAGPA; see 
Section [IV] for further discussion on this procedure) (Fig. [2](a)), by the procedure developed 
in BAPA (Fig. [2tb)), and by our own procedure (Fig. [2l^c)). We have presented in the three 
panels the data pertaining to the Hector Mine earthquake and its aftershocks/recurrences, 
as perceived in the three procedures. One can see that our procedure picks up far more 
events since no time window applies to it. Also, the procedure of DAGPA, since it picks up 
only the record breaking events has minimal links. 
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III. RESULTS 



A. Degree distribution 

The out-degree of a node in these networks describes the number of its recurrences and 
the in-degree describes the number of events which count the selected node as a recurrence 
to themselves. 

Generally speaking, the maximum out-degree (kmax) obtained in the network (Fig. EU^a)) 
is at least an order of magnitude higher than the maximum in-degree (jmax) obtained 
(Fig. El^b)), for all rrith and Cth values. It is clear that large earthquakes will tend to have 
more events associated with them through higher correlation values since Eq. H] is positively 
correlated with magnitude. On the other hand, lower magnitude events require small h/s to 
obtain larger correlation values, i.e. the density of points become important for such events. 
Since jmax has, consistently, lower values compared to k^ax values, this suggests that jmax 
is dependent on lij and, thus, the density of points. Given that the event distribution is 
a power law with respect to magnitude (Eq. [T]) as well as correlation (Fig. [T]), the above 
pattern is sustained with all rrith and Cth values. To summarize, a hub structure is present 
as far as out-degree is concerned while it is absent with the in-degree distribution. 

As Cth is increased and the network becomes sparser, jmax and kmax fall off as power laws, 
for all rrith values, jmax falls off faster and, at Cth = 10^°, the difference with the out-degree 
value is about two orders of magnitude. This again confirms the density dependence of 
the in-degree distribution. Furthermore, at this higher Cth value, jmax becomes magnitude 
independent as well, and all the graphs for the different values converge to the same 
jmax- Similar plots for kmax remain parallel, while falling off, for the different mth values. 
This indicates that the hub structure survives even when the network becomes sparser. 
This is understandable since larger magnitude events dominate kmax distribution and they 
populate the sparse networks produced by the higher cth values. Sparsity induced by higher 
rrith or Cth thresholds does affect the rate of fall-off of kmax- There is a region with a a slow 
fall off (~ —0.17) in the beginning followed by an accelerated fall off (~ —0.65), induced by 
sparsity. The latter sets in at Cth ~ 10^ (Fig. [3](a)). 

The general out-degree distribution, A^^ vs. k (Fig. [3]^c)), where Nk is the number of 
nodes with degree k, mirrors the above scenario seen with maximum out-degree. A power 
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law behavior is obtained for all correlation thresholds and the range of the power law is 
almost over the entire range of out-degree values, i.e. over nearly three decades. The 
exponent of the power law is ~ —2.0. Nevertheless, when k is small, the order of magnitude 
of is practically the same for all correlation thresholds. Since these small k values cannot 
be asssociated with larger magnitude events, this indicates that the lower magnitude events 
have their k values determined by the density of events in the clusters. 

In the case of the in-degree distribution (Fig. Et^d)), the range of degree strengths is much 
smaller and gets even smaller with an increase in cth- This again confirms that sparsity 
affects in-degree significantly since it depends mostly on the density of events. A significant 
plateau is obtained in the middle of the graph, corresponding to intermediate j values, for 
most Cth values, which can be associated with the smaller magnitude events that abound 
in the network. This indicates that a uniform distribution of j values, prescribed by the 
density of events, is the dominant characteristic of the in-degree distribution. 



B. Distribution of recurrence times 



An important feature of earthquake data is that aftershocks decay according to the Omori 



law 



23|,l32| 



n{t) 



[c + ty 

where c and K are constants in time but depend on magnitude m and, p ^ 1. is the 
number of aftershocks per unit time. Earthquakes of all magnitudes have aftershocks which 
decay according to the Omori law. It is widely recorded in literature that even intermediate 
magnitude events can have aftershocks that persist upto years [l^. We investigate this 
by studying the distribution of recurrence times (also referred to as 'waiting times' in the 
literature [a, ll^]) in the data. We define recurrence time here as the time between successive 
events in a cluster surrounding an event. The recurrence times in all the clusters ranges 
from seconds to tens of years. A log-log plot of the distribution is given in Fig. H] where we 
have left out recurrence times < 180 seconds because of the likely error margins. The plot 
shows that a power law holds for most of the range of recurrence times and, in particular, 
for about 5-6 decades. Furthermore, the law holds equally well for all ntth values as the 
overlap of the graphs for mth = 2.5 and ntth = 4.0 shows. The higher recurrence times 
do show some scatter, but that is understandable given the sparsity of data in that region. 
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FIG. 3: The linear variation in the log-log plot of maximum out-degree (kmax) with cth (Fig. [3)[a)) 
shows that k^ax decreases as a power law as Cfh increases, with an accelerated fall off after q/j 
= 10^. The different graphs are for different mth values and it is seen that the behavior is similar 
with the graphs for different rrith values being parallel. The maximum in-degree (jmax) (Fig. El^b)) 
is seen to be couple of magnitudes lower than kmax (cf. Fig. [3{a)). jmax also decreases with cth 
as a power law, as with kmax- However, it also becomes magnitude independent for larger cth 
values as is seen in the figure where all the graphs for different values converge to a common 
value at cth = 10^*^. The general out degree {k) distribution is similar with a power law observed 
(Fig. mc)) there too; the plots are for a fixed rtith (= 2.5) and for different cth values (shown by the 
different symbols), and the power law exponent is ~ —2.0. In Fig. [3]^d), we see that the in-degree 
(j) distribution, on the other hand, has a plateau for the intermediate j values indicating a uniform 
distribution on an average. 
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There is, however, a small change in slope for the higher recurrence times. The initial slope 
of —1.0 changes over to —1.2 for one or two decades before sparsity of data take over. The 
range over which the power law holds seems to increase with increase in mth as seen from 
the graph for mth = 4.0 when compared to the graph for mth = 2.5. We find it remarkable 
that a power law holds for the recurrence times just as the Omori law for the aftershocks, 
which are a subset of the recurrences. In particular, note that the power law is very good up 
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FIG. 4: Recurrence time distribution for cth = 10^ and for two different rrith values (2.5 and 4.0) 
is shown in log-fog plots. Normalized probabilites have been scaled with respect to bin sizes in the 
plot. Recurrence times < 180 sees have not been plotted. The plots show a power law variation 
with an exponent of about —1 for most of the range on the x axis. The range over which the 
plots have a slope of about —1 increases with increasing mfh. This has been verified with different 
mth values and is also seen here from the two graphs for the different rrith values. However, plot 
also becomes more noisy as rrith is increased due to data paucity. Range of the power law remains 
approximately the same for different cth values. 
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to T ^ 10^ s (equivalently t ^ 4 months). This range in time may be different for different 
seismic regions. Considering that there is no objective criterion for fixing the time window 
for identifying aftershocks, we propose that the upper hmit of the time window for the same 
be fixed at the point where the power law fails in the recurrence time plots such as Fig. |H 



C. Distribution of recurrence lengths 

By recurrence length, we mean the distance from an event to each of the other events in 
the same cluster. In other words, we are measuring the scatter in the way the subsequent 
events in the cluster to an event are distributed around the event in question. While there 
is no established law like the Omori law for recurrence len gth s, studies have indicated that 



a power law seems to hold for recurrence lengths as well 1^. To investigate this, we use 



the great circle distance r^j for computing the distance between two event locations. This 



is calculated using the Haversine formula [28|: if (0^, A^) and (0/,A/) are the (latitude, 
longitude) values for the two event locations and, A0 = (pf — (pg and A A = A/ — A^, then 



= RAa 
10 



FIG. 5: The recurrent length distribution for cth = 10'' is shown in Fig. [5^ and for cth = 10^ is 
shown in Fig. [5)3. Both are log-fog plots. Again, normalized probabilites have been scaled with 
respect to bin sizes in the plot. The different symbols corresspond to different values. All the 
plots are seen to be unimodal and peaks at characteristic r values, denoted as r^, that depend 
only on rrith, being larger for higher values of rrith- The exponent of the power law for the part 
before the peak is independent of rrith but depends on cth and goes from 0.45 for cth = 10^ to 0.55 
for Cth = 10^. The exponent of the power law for the part after the peak is also independent of 
rrith, but again depends on cth', it goes from —1.2 for cth = 10'' to —1.08 for cth = 10^. There is, 
understandably, some scatter for the higher r values because of data paucity. As Cth is increased, 
the peaks shift to lower and lower r values until, for cth = lO^'^, it disappears altogether when it 
falls below 100 m; note that we have left out r values < 100 m due to accuracy constraints. By 
rescaling x and y axes suitably, all the graphs can be made to collapse to a single curve as shown 
in the insets. The rescaling function can be represented as Pm{f) = 10~"™'"'-F(r/10""**'") where 
a = 0.51. a is the same for all Cth values. If the peaks for a particular rrith is plotted against cth 
(not shown here), we find again a power law with an exponent of —0.5 which is the same for all 
rrith values. 



OX) 
© 



•nth = 2-5 + 

nijjj = 3.5 ^ 
m,h = 4.0 □ 



S -1 



12 3 

logr 



OX) 

o 



m,h = 2.5 + 
m,h = 3.0 

nijjj = 3.5 'M 

m,h = 4.0 □ 



f sg 



* + 1 



S -2 



1 

logr 



where Aa, the spherical distance, or central angle, between the two event locations is 



Act = 2 arcsin 



' sm^ ( ^ j + cos(0s) cos(0/) sin^ 



R, the radius of the earth, is taken as 6367 km. 

The probability distribution of recurrence lengths is unimodal in nature and peaks at a 
recurrence length (r^) which is larger for increasing values of rrith, for a fixed Cth (Fig. E])- 
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Note that this is consistent with the usual practice of using larger space-time windows to 
collect aftershocks of larger magnitude events. Keep in mind that the distribution with 
lower rrith value contains the data points appearing in the distribution with higher rrith 
values. Also, not all the data points appearing with the lower rrith value is present with 
the higher rrith value. This points to the presence of a characteristic associated with a 
particular m value and this clearly is larger for larger m values. Since the rupture is larger 

with larger earthquakes, we may think of as a typical rupture length associated with 

I ! 

events of that magnitude (see also [14|)- 

varies with Cth as well and the peaks for a fixed rrith shifts to lower values as Cth 
is increased. As a result, the earlier part of the graph where it rises to the peak gets 
progressively shortened as Cth increases and the network becomes sparser. Eventually, for 
Cth = 10^, we see only the part of the graph after the peak. Note that we are plotting only 
pairs of events with > 100 m. As such, the peaks disappear for larger values of Cth since 
the peak values are now < 100 m. The behaviour seen with is completely in accordance 
with Eq. [3l As Cth increases, we are selecting event pairs of larger magnitudes separated by 
smaller distances. The fact that gets smaller with higher cth suggests also that there is 
a range of rupture lengths associated with any magnitude, and not a very sharply defined 
rupture length. Both the rise to the peak value as well as the subsequent fall off are in the 
form of power laws. The slope of the rise to the peak varies with Cth, but stays same for 
different rrith values for a fixed Cth- 

While the graphs for different rrith are separated before the rise to the peak, they are 
collapsed onto a single curve after their peaks. The latter falls off as a power law with an 
exponent ~ —1.2 before data paucity forces scatter. The exponent reduces to —1.08 as Cth is 
increased and number of events decreases. The plots for the different values of rrith collapse 
on to a single curve when the x axis is rescaled as Z/10°™*'' and the y axis by lO"™**" Pm,(r) 
with a = 0.51. The inset plots in Fig. |5] shows the data collapse after the rescaling. Note 
that the value of a remains the same for all Cth values. The value of the rescaling parameter 
a was arrived at by studying the dependence of with nith- Since both rrith and Cth are 
arbitrarily chosen values, the robustness of a, as the latter is changed, shows that the 
values are reliably related across the scale invariant structure of the network. Our studies 



with other data sets 
and is not universal. 



22| show that the value of a depends on the seismicity of the region 
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IV. DISCUSSION AND CONCLUSIONS 



Network analysis of earthquake events have grown in strength and many different ap 



15|). 



proaches have been attempted in constructing the networks (see, for example, 

m 

Nevertheless, some of the conclusions from these different approaches are the same. For ex- 
ample, all are agreed that the out degree distribution is scale free and has a hub structure. 
However, there are also differences in the conclusions. 

In the earlier work of DAGPA, recurrences were studied using a network approach. How- 
ever, the emphasis was on record breaking events which were defined as those which come 
closer to the present one than any earlier. The event that occurred next in time to the 
present one was taken as the first record and, subsequently, record breaking events were 
used to build up the set of records (equivalently, recurrences). In particular, the magni- 
tude of a subsequent event was not taken into consideration and the recurrences were taken 
strictly on spatial closeness criterion. It is possible to construct a minimal network of recur- 
rences by stipulating that only record breaking events qualify as recurrences. Even though 
the definition of immediate successor in time to the present event as a recurrence is not sat- 
isfactory, nevertheless a minimal network will eventually discard those events when a record 
breaking event is found. 

The data set used in DAGPA is almost the same as the one used in this present study 
and the distribution of recurrence lengths in DAGPA turns out to be similar to the one here 
(Fig. |5]). This is not surprising because it is a minimal network of recurrences, as compared 
to the full set of recurrences considered here. A characteristic length emerges, with the 
minimal network, at which the distribution peaks which was identified in DAGPA to be 
a typical rupture length associated with that magnitude threshold. The slope of the part 
after the peak has the same slope as what is obtained in the present study. However, the 
part before the peak was claimed to have a slope of unity confirming a random distribution 
of events within a disc of 100 m (the location error bar in the data set) around any event. 
However, we find a power law behaviour in this part too with a growth slope ~ 0.5. In 
particular, for each Cth, the graphs for different rrith have the same slope and are parallel. As 
Cth is increased, this part is shortened and eventually disappears for cth = 10^°. Since only 
larger magnitude events are present with these high cth values, this points to the fact that 
the big events do not cluster at closeby distances. On the other hand, this also implies that 
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small magnitude events predominantly cluster at nearby distances. DAGPA reports that 
no systematic errors have been reported in literature for magnitude dependence of location 
errors, and argues that the location errors which are < 100 m, tend to randomize the 
distribution for recurrence lengths < 100 m. However, the order we see in the portion before 
the peak in Fig. El and similar plots for different Cth values, suggests that location errors 
could be magnitude dependent and should be looked for. It is, of course, to be emphasized 
that our analysis includes only the correlated earthquakes which brings in the order seen in 
the portion before the peak. 

It is instructive to compare our results for recurrence times with another study on 'waiting 
times'. Bak et al. 12|] had analyzed data from Southern California and shown that a unified 
scaling law exists for waiting times. They had divided the study area into grids of varying 
sizes and compared the waiting time distribution across the grid. They showed that a scaling 
law of the form 

T"Ps,l{T) = f (TS-^L'^f) 

holds where P is the waiting time distribution, T the waiting time, m = log S and L is 
the linear size of the grid. The function / was seen to consist of a constant part and a 
decaying part, separated by a sharp kink. They concluded that the constant part consisted 
of the correlated earthquakes and the decaying part with the uncorrelated earthquakes, i.e. 
independendent earthquakes seeded by the nonzero driving rate of plate tectonics. Several 



other studies have since appeared ana 
Eq. [2] as well as the 'waiting times 



,^zin g th is unified scaling law which combines Eq. [H 

Our studies are based on the same area and practically the same data set. However, we 
have divided the events into clusters based on correlation values. In other words, we have 
picked only the correlated earthquakes for analysis. There is, hence, no recurrence times 
coming from uncorrelated earthquakes. Therefore, there is no kink in the data as seen in 
Bak et al. 12| and our plot has only one slope which agrees qualitatively with the portion 



of the graph in Bak et al. 



12j with correlated earthquakes. 



As already mentioned, our definition of the correlation, Eq. [3] and Eq. HJ differs from 
BAPA in that we have not included tij in our definition and, also, we evaluate the correla- 
tions going forwards in time as against the backward associations in BAPA. The approach 
developed in BAPA is able to address the question of identifying aftershocks in an unbiased 



14 



manner and without imposing pre-determined space-time windows. On the other hand, we 
expect the present approach, coupled with the induction of tij, to identify foreshocks, which 
we intend to pursue in a future work. Finally, we are also proposing that the upper limit 
of the time window where the power law holds in recurrence time distribution plots such as 
Fig. m could be used as an objective criterion to fix the limit of the aftershock window. 
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